% =========================================================================
% tmp.m
% Function for:
%
% STICKY PRICES AND MONETARY POLICY: EVIDENCE FROM DISAGGREGATED U.S. DATA
% by Jean Boivin, Marc Giannoni, and Ilian Mihov 
% American Economic Review
% 
% Called by irfbootfac.m
% Computes PC factors
% Syntax:
%   [factors, lam, ebp]= tmp(datam,dataq,nfac,fm)
% 
% where the inputs are:     
%   datam = Xt vector in the observation equation
%   dataq = matrix of quarterly series
%   nfac = number of factors
%
% and the outputs are:
%   factors = estimated PC factors
%   lam = estimated PC factors coeficients
%   ebp = estimated error terms
% =========================================================================

function [factors, lam, ebp]= tmp(datam,dataq,nfac,fm);

% Initial estimates of factors from monthly data

[t nb]=size(datam);
 xx=datam'*datam;
[evec,eval]=eig(xx);

% Sorting evec so that they correspond to eval in descending order
[eval,index]=sort(diag(eval));      % sorts eigenvalues in an ascending order 
index=flipud(index); 		        % to get descending order
evc=zeros(nb,nb);                   
for i=1:nb                          
    evc(:,i)=evec(:,index(i));
end

% Estimate the factors, coeficients and error terms
lam = sqrt(nb)*evc(:,1:nfac);		% lam is matrix nb*nfac
factors=datam*lam/nb;				% fac is t*nfac
ebp=datam-(factors*lam');
ssr=sum(sum(ebp.^2));

%Adjustments in the case of quarterly series

if size(dataq,1) > 1
    % Recursions on complete data set using factors from monthly data as starting value  
    [a,invaa] = matrixa(t,fm);
    for it=1:60   
        % Step 1: Construct fitted values based on most recent estimate of the factors  
        xq=repmat(NaN,[size(datam,1) size(dataq,2)]);
        qy=[];
        ssrq=0;
        for j=1:size(dataq,2);
            [dataqj,qyj,iq,ssrqj] = estepmf(dataq(:,j),factors,a,invaa,21,fm);
            ssrq=ssrq+ssrqj;
            xq(:,j)=dataqj;
            qy=[qy qyj];
        end
        
        % Step 2: Update factors using fitted values of the mixed frequencies panel
        xtot=[datam xq];
        xxtot=xtot'*xtot;
        [evec,eval]=eig(xxtot);
        [eval,index]=sort(diag(eval));
        index=flipud(index);
        evc=zeros(size(evec));
        for i=1:size(evec,1)
            evc(:,i)=evec(:,index(i));
        end
        lam = sqrt(size(evec,1))*evc(:,1:nfac);
        factors=xtot*lam/size(evec,1);
        ebp=xtot-(factors*lam');
        
        % observe convergence with the sum of squared residuals for not missing observations, (~A) has 0's at missing obs
        v=(ssrq)/size(xtot,2);     
    end
end